function[H]=CM(u1,u2,v,theta,D,k,xi,N,B)
%Continuum Model Hamiltonian
%basis {1A,1B,2A,2B}

% k=[kx;ky];

N=2*N;

w=exp(1i*2*pi/3);

G1=8*pi/sqrt(3)*sin(theta/2)*[-1/2;-sqrt(3)/2];
G2=8*pi/sqrt(3)*sin(theta/2)*[1/2;-sqrt(3)/2];

%Interlayer coupling
U0=[u1,u2;
    u2,u1];% zero momentum transfer

U1=[u1,u2*w^(-xi);
    u2*w^(xi),u1];% momentum transfer G1

U2=[u1,u2*w^(xi);
    u2*w^(-xi),u1];% momentum transfer G2

H=kron(speye(N),U0)+kron(diag(ones(N-1,1),1),U1);
% H=H+H';
H=kron(speye(N),H)+kron(diag(ones(N-1,1),1),kron(speye(N),U2));  
% H=H+H';
%Intralyer Dirac Hamiltonian
H1=0*H;
H2=0*H;
for n=1:N
    for m=1:N
S=zeros(N^2);
S(m+(n-1)*N,m+(n-1)*N)=1;
kk=k-xi*N*(G1+G2)/2+(m-1)*xi*G1+(n-1)*xi*G2;
[h1,h2]=DH(v,theta,D,kk,xi,B);
H1=H1+kron(S,h1);
H2=H2+kron(S,h2);
    end
end


H=[H1,H;
   H',H2];

H=full(H);

end